home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / tutorials / geometer / constraint.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  16.0 KB  |  556 lines
  1. /*
  2.  * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17.  
  18. #include "parse.h"
  19. #include <math.h>
  20. #include <stream.h>
  21. #include <sys/types.h>
  22. #include <sys/stat.h>
  23. #include "generic.h"
  24. #include "gizmo.h"
  25.  
  26. // vertonlineline(): makes the vertex v at the intersection of
  27. // lines l1 and l2.
  28.  
  29. void vertonlineline(vertex *v, line *l1, line *l2)
  30. {
  31.     v->xw = l1->YW*l2->W - l1->W*l2->YW;
  32.     v->yw = l1->W*l2->XW - l1->XW*l2->W;
  33.     v->w = l1->XW*l2->YW - l1->YW*l2->XW;
  34.     if (v->w != 0.0) {
  35.     v->xw /= v->w; v->yw /= v->w; v->w = 1.0;
  36.     }
  37. }
  38.  
  39. float vvdist(vertex *v1, vertex *v2)
  40. {
  41.     if (v1->w == 0.0 || v2->w == 0) return 1.0e30;
  42.     float x1 = v1->xw;
  43.     float x2 = v2->xw;
  44.     float y1 = v1->yw;
  45.     float y2 = v2->yw;
  46.     return sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
  47. }
  48.  
  49. float vldist(vertex *v, line *l)    // perp distance: v to l
  50. {
  51.     if (v->w == 0 || (l->XW == 0.0 && l->YW == 0.0)) return 1.0e30;
  52.     return fabs((l->XW*v->xw+l->YW*v->yw+l->W)/sqrt(l->XW*l->XW+l->YW*l->YW));
  53. }
  54.  
  55. void solvevlc(line *l, circle *c, vertex *v, long rootnumber)
  56. {
  57.     float disc, vx, vy, wx, wy, alpha;
  58.     long ok = 0;
  59.  
  60.     if (l->v1.w != 0.0) {
  61.     vx = l->v1.xw;
  62.     vy = l->v1.yw;
  63.     ok = 1;
  64.     } else {
  65.     vx = l->v1.xw + l->v2.xw;
  66.     vy = l->v1.yw + l->v2.yw;
  67.     }
  68.     if (l->v2.w != 0.0) {
  69.     wx = l->v2.xw;
  70.     wy = l->v2.yw;
  71.     ok = 1;
  72.     } else {
  73.     wx = l->v1.xw + l->v2.xw;
  74.     wy = l->v1.yw + l->v2.yw;
  75.     }
  76.     float cx = c->center.xw;
  77.     float cy = c->center.yw;
  78.     float R = c->radius;
  79.     float A = (vx-wx)*(vx-wx)+(vy-wy)*(vy-wy);
  80.     float B = 2*((vx-wx)*(wx-cx)+(vy-wy)*(wy-cy));
  81.     float C = (wx-cx)*(wx-cx)+(wy-cy)*(wy-cy) - R*R;
  82.     if ((disc = B*B-4.0*A*C) < 0.0 || (ok == 0)) {
  83.     v->xw = v->yw = v->w = 0.0;
  84.     return;
  85.     }
  86.     if (rootnumber == 1) alpha = ((-B + sqrt(disc))/(2.0*A));
  87.     else alpha = ((-B - sqrt(disc))/(2.0*A));
  88.     v->xw = alpha*(vx)+(1-alpha)*(wx);
  89.     v->yw = alpha*(vy)+(1-alpha)*(wy);
  90.     v->w = 1.0;
  91. }
  92.  
  93. // homogenizeline(): assumes that the two vertices are correct;
  94. // figures out the line's homogeneous coordinates from those.
  95.  
  96. void homogenizeline(line *l)
  97. {
  98.     if ((l->v1.xw == 0.0 && l->v1.yw == 0.0 && l->v1.w == 0.0) ||
  99.         (l->v2.xw == 0.0 && l->v2.yw == 0.0 && l->v2.w == 0.0)) {
  100.         l->XW = l->YW = l->W = 0.0;
  101.         return;
  102.     }
  103.     l->XW = l->v1.yw*l->v2.w - l->v1.w*l->v2.yw;
  104.     l->YW = l->v2.xw*l->v1.w - l->v1.xw*l->v2.w;
  105.     l->W = l->v1.xw*l->v2.yw - l->v2.xw*l->v1.yw;
  106. }
  107.  
  108. void line_vertcircletan(line *l)
  109. {
  110.     l->v1.w = l->v2.w = 1.0;
  111.     vertex *v = (vertex *)l->c.p1;
  112.     circle *c = (circle *)l->c.p2;
  113.     if ((c->c.type == CIRCVERTVERT && c->c.p2 == v) ||
  114.         (v->c.type == VERTONCIRCLE && v->c.p1 == c) ||
  115.     (v->c.type == VERTLINECIRC1 && v->c.p2 == c) || 
  116.     (v->c.type == VERTLINECIRC2 && v->c.p2 == c) || 
  117.     (v->c.type == VERTCIRCCIRC1 && (v->c.p1 == c || v->c.p2 == c)) || 
  118.     (v->c.type == VERTCIRCCIRC2 && (v->c.p1 == c || v->c.p2 == c)) ||
  119.     (c->c.type == CIRCVERTVERTVERT && ((c->c.p1 == v) ||
  120.         (c->c.p2 == v) || (c->c.p3 == v)))) {
  121.         // Vertex is defined to be on the circle.
  122.         l->v1.xw = v->xw;
  123.         l->v1.yw = v->yw;
  124.         if (v->yw == c->center.yw) {
  125.         l->v2.xw = v->xw;
  126.         l->v2.yw = 0.0;
  127.         } else {
  128.             l->v2.xw = v->xw + 1.0;
  129.         l->v2.yw = v->yw - (v->xw - c->center.xw)/(v->yw - c->center.yw);
  130.         }
  131.         homogenizeline(l);
  132.         return;
  133.     }
  134.     // XXX assumes finite center
  135.     float qx = v->xw - c->center.xw;
  136.     float qy = v->yw - c->center.yw;
  137.     float R = sqrt(qx*qx+qy*qy);
  138.     float r = c->radius;
  139.     float xx = r*r/R;
  140.     float d = 1.0-xx/R;
  141.     if (d < 0.0) {
  142.     l->v1.xw = l->v2.xw = 0.0;
  143.     l->v1.yw = 0.0; l->v2.yw = 0.0;
  144.     l->v1.w = l->v2.w = 0.0;
  145.     homogenizeline(l);
  146.     return;
  147.     }
  148.     float yy = r*sqrt(d);
  149.     if (l->c.type == LINEVERTCIRC2) yy = -yy;
  150.     float xx1 = xx*qx/R - yy*qy/R;
  151.     float yy1 = xx*qy/R + yy*qx/R;
  152.     l->v1.xw = v->xw; l->v1.yw = v->yw; l->v1.w = v->w;
  153.     l->v2.xw = xx1 + c->center.xw;
  154.     l->v2.yw = yy1 + c->center.yw;
  155.     homogenizeline(l);
  156. }
  157.  
  158. void line_circcircleint(line *l)
  159. {
  160.     // XXX assumes circle's centers are finite
  161.     l->v1.w = l->v2.w = 1.0;
  162.     circle *c1 = (circle *)l->c.p1;
  163.     circle *c2 = (circle *)l->c.p2;
  164.     float c1x = c1->center.xw, c1y = c1->center.yw;
  165.     float c2x = c2->center.xw, c2y = c2->center.yw;
  166.     float r1 = c1->radius, r2 = c2->radius;
  167.     if (r1 == r2) return;
  168.     float vx = c2x - (c2x-c1x)*(r2/(r1+r2));
  169.     float vy = c2y - (c2y-c1y)*(r2/(r1+r2));
  170.     float qx = vx - c1->center.xw;
  171.     float qy = vy - c1->center.yw;
  172.     float R = sqrt(qx*qx+qy*qy);
  173.     float r = c1->radius;
  174.     float xx = r*r/R;
  175.     float d = 1.0-xx/R;
  176.     float qx2 = vx - c2->center.xw;
  177.     float qy2 = vy - c2->center.yw;
  178.     float R2 = sqrt(qx2*qx2+qy2*qy2);
  179.     float rr2 = c2->radius;
  180.     float xx2 = rr2*rr2/R2;
  181.     float d2 = 1.0-xx2/R2;
  182.     if (d < 0.0) {
  183.     l->v1.xw = l->v2.xw = 0.0;
  184.     l->v1.yw = 0.0; l->v2.yw = 0.0;
  185.     l->v1.w = l->v2.w = 0.0;
  186.     homogenizeline(l);
  187.     return;
  188.     }
  189.     float yy = r*sqrt(d);
  190.     if (l->c.type == LINECIRCCIRCINT2) yy = -yy;
  191.     float yy2 = rr2*sqrt(d2);
  192.     if (l->c.type == LINECIRCCIRCINT2) yy2 = -yy2;
  193.     float xx1 = xx*qx/R - yy*qy/R;
  194.     float yy1 = xx*qy/R + yy*qx/R;
  195.     float xxx2 = xx2*qx2/R2 - yy2*qy2/R2;
  196.     float yyy2 = xx2*qy2/R2 + yy2*qx2/R2;
  197.     //l->v1.xw = vx; l->v1.yw = vy;
  198.     l->v2.xw = xx1 + c1->center.xw;
  199.     l->v2.yw = yy1 + c1->center.yw;
  200.     l->v1.xw = xxx2 + c2->center.xw;
  201.     l->v1.yw = yyy2 + c2->center.yw;
  202.     homogenizeline(l);
  203. }
  204.  
  205. void line_circcircleext(line *l)
  206. {
  207.     // assumes finite circle centers
  208.     l->v1.w = l->v2.w = 1.0;
  209.     circle *c1 = (circle *)l->c.p1;
  210.     circle *c2 = (circle *)l->c.p2;
  211.     float c1x = c1->center.xw, c1y = c1->center.yw;
  212.     float c2x = c2->center.xw, c2y = c2->center.yw;
  213.     float r1 = c1->radius, r2 = c2->radius;
  214.     if (r1 == r2) return;
  215.     float vx = c2x + (c2x-c1x)*(r2/(r1-r2));
  216.     float vy = c2y + (c2y-c1y)*(r2/(r1-r2));
  217.     float qx = vx - c1->center.xw;
  218.     float qy = vy - c1->center.yw;
  219.     float R = sqrt(qx*qx+qy*qy);
  220.     float r = c1->radius;
  221.     float xx = r*r/R;
  222.     float d = 1.0-xx/R;
  223.     float qx2 = vx - c2->center.xw;
  224.     float qy2 = vy - c2->center.yw;
  225.     float R2 = sqrt(qx2*qx2+qy2*qy2);
  226.     float rr2 = c2->radius;
  227.     float xx2 = rr2*rr2/R2;
  228.     float d2 = 1.0-xx2/R2;
  229.     if (d < 0.0) {
  230.     l->v1.xw = l->v2.xw = 0.0;
  231.     l->v1.yw = 0.0; l->v2.yw = 0.0;
  232.     l->v1.w = l->v2.w = 0.0;
  233.     homogenizeline(l);
  234.     return;
  235.     }
  236.     float yy = r*sqrt(d);
  237.     if (l->c.type == LINECIRCCIRCEXT2) yy = -yy;
  238.     float yy2 = rr2*sqrt(d2);
  239.     if (l->c.type == LINECIRCCIRCEXT2) yy2 = -yy2;
  240.     float xx1 = xx*qx/R - yy*qy/R;
  241.     float yy1 = xx*qy/R + yy*qx/R;
  242.     float xxx2 = xx2*qx2/R2 - yy2*qy2/R2;
  243.     float yyy2 = xx2*qy2/R2 + yy2*qx2/R2;
  244.     //l->v1.xw = vx; l->v1.yw = vy;
  245.     l->v2.xw = xx1 + c1->center.xw;
  246.     l->v2.yw = yy1 + c1->center.yw;
  247.     l->v1.xw = xxx2 + c2->center.xw;
  248.     l->v1.yw = yyy2 + c2->center.yw;
  249.     homogenizeline(l);
  250. }
  251.  
  252. void vert_circcircle(vertex *v)
  253. {
  254.     circle *c1 = (circle *)v->c.p1;
  255.     circle *c2 = (circle *)v->c.p2;
  256.     float c1x = c1->center.xw, c1y = c1->center.yw;
  257.     float c2x = c2->center.xw, c2y = c2->center.yw;
  258.     float r1 = c1->radius, r2 = c2->radius;
  259.     float Q = sqrt((c1x-c2x)*(c1x-c2x)+(c1y-c2y)*(c1y-c2y));
  260.     float Rcost = -(r2*r2 - r1*r1 - Q*Q)/(2*Q);
  261.     float d = r1*r1 - Rcost*Rcost;
  262.     v->w = 1.0;
  263.     if (d < 0.0) {
  264.     v->xw = v->yw = v->w = 0.0;
  265.     return;
  266.     }
  267.     float Rsint = sqrt(d);
  268.     float vx = c1x + Rcost*(c2x-c1x)/Q;
  269.     float vy = c1y + Rcost*(c2y-c1y)/Q;
  270.     if (v->c.type == VERTCIRCCIRC1) {
  271.     v->xw = vx + Rsint*(c2y-c1y)/Q;
  272.     v->yw = vy + Rsint*(c1x-c2x)/Q;
  273.     } else {
  274.     v->xw = vx - Rsint*(c2y-c1y)/Q;
  275.     v->yw = vy - Rsint*(c1x-c2x)/Q;
  276.     }
  277. }
  278.  
  279. void vert_lineconic(vertex *v)
  280. {
  281.     line *l = (line *)v->c.p1;
  282.     conic *c = (conic *)v->c.p2;
  283.     float A = l->XW, B = l->YW, C = l->W;
  284.     if (A != 0.0) {
  285.     float AA = c->aa*B*B/(A*A) - c->bb*B/A + c->cc;
  286.     float BB = 2.0*c->aa*B*C/(A*A) - c->bb*C/A - c->dd*B/A + c->ee;
  287.     float CC = c->aa*C*C/(A*A) - c->dd*C/A + c->ff;
  288.     float D = BB*BB - 4.0*AA*CC;
  289.     if (D < 0) {
  290.         v->xw = 0.0; v->yw = 0.0; v->w = 0.0;
  291.         return;
  292.     }
  293.     if (v->c.type == VERTLINECONIC1) {
  294.         v->yw = (-BB + sqrt(D))/(2.0*AA);
  295.         v->xw = -B*v->yw/A - C/A;
  296.         v->w = 1.0;
  297.     } else {
  298.         v->yw = (-BB - sqrt(D))/(2.0*AA);
  299.         v->xw = -B*v->yw/A - C/A;
  300.         v->w = 1.0;
  301.     }
  302.     } else {
  303.     float AA = c->aa;
  304.     float BB = c->dd - C*c->bb/B;
  305.     float CC = c->cc*C*C/(B*B) - c->ee*C/B + c->ff;
  306.     float D = BB*BB - 4.0*AA*CC;
  307.     if (D < 0) {
  308.         v->xw = 0.0; v->yw = 0.0; v->w = 0.0;
  309.         return;
  310.     }
  311.     if (v->c.type == VERTLINECONIC1) {
  312.         v->yw = -C/B;
  313.         v->xw = (-BB + sqrt(D))/(2.0*AA);
  314.         v->w = 1.0;
  315.     } else {
  316.         v->yw = -C/B;
  317.         v->xw = (-BB - sqrt(D))/(2.0*AA);
  318.         v->w = 1.0;
  319.     }
  320.     }
  321. }
  322.  
  323. void circle_lll(circle *c)
  324. {
  325.     vertex v1, v2, v3;
  326.     vertex w1, w2;
  327.     long circindex;
  328.     line ll1, ll2;
  329.     
  330.     line *l1 = (line *)c->c.p1;
  331.     line *l2 = (line *)c->c.p2;
  332.     line *l3 = (line *)c->c.p3;
  333.     vertonlineline(&v3, l1, l2);
  334.     vertonlineline(&v2, l1, l3);
  335.     vertonlineline(&v1, l2, l3);
  336.     switch (c->c.type) {
  337.     case CIRCLINELINELINE1: circindex = 0; break;
  338.     case CIRCLINELINELINE2: circindex = 1; break;
  339.     case CIRCLINELINELINE3: circindex = 2; break;
  340.     case CIRCLINELINELINE4: circindex = 3; break;
  341.     }
  342.     float R12 = vvdist(&v1, &v2);
  343.     float R13 = vvdist(&v1, &v3);
  344.     float R23 = vvdist(&v2, &v3);
  345.     if (R12 == 0.0 || R13 == 0.0 || R23 == 0.0) {
  346.     c->center.xw = c->center.yw = 1.0e30;
  347.     c->radius = 1.0;
  348.     return;
  349.     }
  350.     float alpha1 = 1 + ((circindex&1) ? R12/R13 : - R12/R13);
  351.     float alpha2 = 1 + ((circindex&2) ? R12/R23 : - R12/R23);
  352.     w1.xw = (alpha1*v1.xw + (1-alpha1)*v3.xw + v2.xw)/2.0;
  353.     w1.yw = (alpha1*v1.yw + (1-alpha1)*v3.yw + v2.yw)/2.0;
  354.     w1.w = 1.0;
  355.     // v1-w1 is one of two angle bisectors
  356.     w2.xw = (alpha2*v2.xw + (1-alpha2)*v3.xw + v1.xw)/2.0;
  357.     w2.yw = (alpha2*v2.yw + (1-alpha2)*v3.yw + v1.yw)/2.0;
  358.     w2.w = 1.0;
  359.     ll1.v1 = v1; ll1.v2 = w1;
  360.     ll2.v1 = v2; ll2.v2 = w2;
  361.     homogenizeline(&ll1); homogenizeline(&ll2);
  362.     vertonlineline(&c->center, &ll1, &ll2);
  363.     c->radius = vldist(&c->center, l1);
  364. }
  365.  
  366. void circle_vvv(circle *c)
  367. {
  368.     c->center.w = 1.0;
  369.     vertex *v1 = (vertex *)c->c.p1;
  370.     vertex *v2 = (vertex *)c->c.p2;
  371.     vertex *v3 = (vertex *)c->c.p3;
  372.     float bx = v1->xw; float by = v1->yw;
  373.     float cx = v2->xw; float cy = v2->yw;
  374.     float dx = v3->xw; float dy = v3->yw;
  375.     float temp = cx*cx+cy*cy;
  376.     float bc = (bx*bx + by*by - temp)/2.0;
  377.     float cd = (temp - dx*dx - dy*dy)/2.0;
  378.     float det = (bx-cx)*(cy-dy)-(cx-dx)*(by-cy);
  379.     if (fabs(det) < 1.0e-6) {
  380.     c->center.xw = c->center.yw = 1.0;
  381.     c->center.w = 0.0;
  382.     c->v1 = *v1;
  383.     c->v2 = *v2;
  384.     c->v3 = *v3;
  385.     return;
  386.     }
  387.     det = 1/det;
  388.     c->center.xw = (bc*(cy-dy)-cd*(by-cy))*det;
  389.     c->center.yw = ((bx-cx)*cd-(cx-dx)*bc)*det;
  390.     cx = c->center.xw; cy = c->center.yw;
  391.     c->radius = sqrt((cx-bx)*(cx-bx)+(cy-by)*(cy-by));
  392. }
  393.  
  394. void circle_ccinvert(circle *c)
  395. {
  396.     c->center.w = 1.0;
  397.     circle *inv = (circle *)c->c.p2;
  398.     circle *crc = (circle *)c->c.p1;
  399.     vertex ctr = inv->center;
  400.     float r = inv->radius;
  401.     float px, py, qx, qy, rx, ry;
  402.     px = qx = crc->center.xw;
  403.     rx = px + crc->radius;
  404.     ry = crc->center.yw;
  405.     py = crc->center.yw + crc->radius;
  406.     qy = crc->center.yw - crc->radius;
  407.     float bx, by, cx, cy, dx, dy;
  408.     float f;
  409.     f = r*r/((px-ctr.xw)*(px-ctr.xw)+(py-ctr.yw)*(py-ctr.yw));
  410.     bx = ctr.xw + (px - ctr.xw) * f;
  411.     by = ctr.yw + (py - ctr.yw) * f;
  412.     f = r*r/((qx-ctr.xw)*(qx-ctr.xw)+(qy-ctr.yw)*(qy-ctr.yw));
  413.     cx = ctr.xw + (qx - ctr.xw) * f;
  414.     cy = ctr.yw + (qy - ctr.yw) * f;
  415.     f = r*r/((rx-ctr.xw)*(rx-ctr.xw)+(ry-ctr.yw)*(ry-ctr.yw));
  416.     dx = ctr.xw + (rx - ctr.xw) * f;
  417.     dy = ctr.yw + (ry - ctr.yw) * f;
  418.     float bc = (bx*bx + by*by - cx*cx - cy*cy)/2.0;
  419.     float cd = (cx*cx + cy*cy - dx*dx - dy*dy)/2.0;
  420.     float det = (bx-cx)*(cy-dy)-(cx-dx)*(by-cy);
  421.     if (det == 0) {
  422.     c->center.xw = c->center.yw = 1.0e30;
  423.     c->radius = 1.0;
  424.     return;
  425.     }
  426.     c->center.xw = (bc*(cy-dy)-cd*(by-cy))/det;
  427.     c->center.yw = ((bx-cx)*cd-(cx-dx)*bc)/det;
  428.     cx = c->center.xw; cy = c->center.yw;
  429.     c->radius = sqrt((cx-bx)*(cx-bx)+(cy-by)*(cy-by));
  430. }
  431.  
  432.  
  433. void v_vcinvert(vertex *v)
  434. {
  435.     circle *inv = (circle *)v->c.p2;
  436.     vertex *orig = (vertex *)v->c.p1;
  437.     vertex ctr = inv->center;
  438.     float r = inv->radius;
  439.     float px, py;
  440.     px = orig->xw;
  441.     py = orig->yw;
  442.     float f;
  443.     f = r*r/((px-ctr.xw)*(px-ctr.xw)+(py-ctr.yw)*(py-ctr.yw));
  444.     v->xw = ctr.xw + (px - ctr.xw) * f;
  445.     v->yw = ctr.yw + (py - ctr.yw) * f;
  446.     v->w = 1.0;
  447. }
  448.  
  449. void circle_linecinvert(circle *c)
  450. {
  451.     c->center.w = 1.0;
  452.     circle *inv = (circle *)c->c.p2;
  453.     line *l = (line *)c->c.p1;
  454.     vertex ctr = inv->center;
  455.     float r = inv->radius;
  456.     float px, py, qx, qy, rx, ry;
  457.     px = l->v1.xw; py = l->v1.yw;
  458.     qx = l->v2.xw; qy = l->v2.yw;
  459.     rx = px + 2*(qx-px);
  460.     ry = py + 2*(qy-py);
  461.     float bx, by, cx, cy, dx, dy;
  462.     float f;
  463.     f = r*r/((px-ctr.xw)*(px-ctr.xw)+(py-ctr.yw)*(py-ctr.yw));
  464.     bx = ctr.xw + (px - ctr.xw) * f;
  465.     by = ctr.yw + (py - ctr.yw) * f;
  466.     f = r*r/((qx-ctr.xw)*(qx-ctr.xw)+(qy-ctr.yw)*(qy-ctr.yw));
  467.     cx = ctr.xw + (qx - ctr.xw) * f;
  468.     cy = ctr.yw + (qy - ctr.yw) * f;
  469.     f = r*r/((rx-ctr.xw)*(rx-ctr.xw)+(ry-ctr.yw)*(ry-ctr.yw));
  470.     dx = ctr.xw + (rx - ctr.xw) * f;
  471.     dy = ctr.yw + (ry - ctr.yw) * f;
  472.     float bc = (bx*bx + by*by - cx*cx - cy*cy)/2.0;
  473.     float cd = (cx*cx + cy*cy - dx*dx - dy*dy)/2.0;
  474.     float det = (bx-cx)*(cy-dy)-(cx-dx)*(by-cy);
  475.     if (det == 0) {
  476.     c->center.xw = c->center.yw = 1.0e30;
  477.     c->radius = 1.0;
  478.     return;
  479.     }
  480.     c->center.xw = (bc*(cy-dy)-cd*(by-cy))/det;
  481.     c->center.yw = ((bx-cx)*cd-(cx-dx)*bc)/det;
  482.     cx = c->center.xw; cy = c->center.yw;
  483.     c->radius = sqrt((cx-bx)*(cx-bx)+(cy-by)*(cy-by));
  484. }
  485.  
  486. void ratio_vvv(length *r)
  487. {
  488.     vertex *v1 = (vertex *)r->c.p1;
  489.     vertex *v2 = (vertex *)r->c.p2;
  490.     vertex *v3 = (vertex *)r->c.p3;
  491.     float dist1 = vvdist(v1, v2);
  492.     float dist2 = vvdist(v1, v3);
  493.     float dot = (v3->xw-v1->xw)*(v2->xw-v1->xw)
  494.                         + (v3->yw-v1->yw)*(v2->yw-v1->yw);
  495.     if (dist2 == 0) r->value = 1.0e30;
  496.     else r->value = dist1/dist2;
  497.     if (dot < 0.0) r->value *= -1.0;
  498. }
  499.  
  500. void vlvmirror(vertex *v)
  501. {
  502.     line *l = (line *)v->c.p1;
  503.     vertex *vv = (vertex *)v->c.p2;
  504.     vertex *v1 = &l->v1;
  505.     vertex *v2 = &l->v2;
  506.     float x0 = v1->xw;
  507.     float y0 = v1->yw;
  508.     float x1 = v2->xw;
  509.     float y1 = v2->yw;
  510.     float px = vv->xw,  py = vv->yw;
  511.     float xw, yw;
  512.     float alpha = ((x1-px)*(x1-x0) + (y1-py)*(y1-y0))/
  513.             ((x1-x0)*(x1-x0)+(y1-y0)*(y1-y0));
  514.     xw = alpha*x0 + (1-alpha)*x1;
  515.     yw = alpha*y0 + (1-alpha)*y1;
  516.     v->xw = 2*xw - px;
  517.     v->yw = 2*yw - py;
  518.     v->w = 1.0;
  519. }
  520.  
  521. void v_vvratio(vertex *v)
  522. {
  523.     vertex *v1 = (vertex *)v->c.p1;
  524.     vertex *v2 = (vertex *)v->c.p2;
  525.     float rat = ((length *)(v->c.p3))->value;
  526.     v->xw = (1-rat)*v1->xw + rat*v2->xw;
  527.     v->yw = (1-rat)*v1->yw + rat*v2->yw;
  528.     v->w = 1.0;
  529. }
  530.  
  531. void a_vvv(length *a)
  532. {
  533.     vertex *v1 = (vertex *)a->c.p1;
  534.     vertex *v2 = (vertex *)a->c.p2;
  535.     vertex *v3 = (vertex *)a->c.p3;
  536.     float x2 = v3->xw - v2->xw;
  537.     float y2 = v3->yw - v2->yw;
  538.     float x1 = v1->xw - v2->xw;
  539.     float y1 = v1->yw - v2->yw;
  540.     a->value = atan2(y2, x2) - atan2(y1, x1);
  541.     if (a->value < 0.0) a->value += 2.0*PI;
  542. }
  543.  
  544. void v_avv(vertex *v)
  545. {
  546.     length *a = (length *)v->c.p1;
  547.     vertex *v1 = (vertex *)v->c.p2;
  548.     vertex *v2 = (vertex *)v->c.p3;
  549.     float x1 = v1->xw - v2->xw;
  550.     float y1 = v1->yw - v2->yw;
  551.     float c = cos(a->value), s = sin(a->value);
  552.     v->xw = v2->xw + x1*c - y1*s;
  553.     v->yw = v2->yw + x1*s + y1*c;
  554.     v->w = 1.0;
  555. }
  556.